Manipulating local coordination of copper single atom catalyst enables efficient CO2-to-CH4 conversion

Electrochemical CO2 conversion to methane, powered by intermittent renewable electricity, provides an entrancing opportunity to both store renewable electric energy and utilize emitted CO2. Copper-based single atom catalysts are promising candidates to restrain C-C coupling, suggesting feasibility in further protonation of CO* to CHO* for methane production. In theoretical studies herein, we find that introducing boron atoms into the first coordination layer of Cu-N4 motif facilitates the binding of CO* and CHO* intermediates, which favors the generation of methane. Accordingly, we employ a co-doping strategy to fabricate B-doped Cu-Nx atomic configuration (Cu-NxBy), where Cu-N2B2 is resolved to be the dominant site. Compared with Cu-N4 motifs, as-synthesized B-doped Cu-Nx structure exhibits a superior performance towards methane production, showing a peak methane Faradaic efficiency of 73% at −1.46 V vs. RHE and a maximum methane partial current density of −462 mA cm−2 at −1.94 V vs. RHE. Extensional calculations utilizing two-dimensional reaction phase diagram analysis together with barrier calculation help to gain more insights into the reaction mechanism of Cu-N2B2 coordination structure.

natural gas, possesses a good compatibility with the existing infrastructure for storage, distribution, and consumption [18][19][20] . Then, electrocatalytic conversion of CO 2 into CH 4 offers an entrancing opportunity to both storing renewable electric energy and utilizing CO 2 emissions. In general, starting with the CO* intermediates of CO 2 RR, protonation of CO* to CHO* leads to CH 4 , whereas the competitive CO* dimerization generates C 2 products [21][22][23][24] . On the other hand, the non-optimized binding of CO* intermediates will result in release of gaseous CO, further suppressing the CH 4 selectivty 25 . To perform CO 2 -to-CH 4 economically at scale, a catalyst capable of mediating the efficient formation of CH 4 with high selectivity and productivity is a prerequisite. Copper (Cu) has been extensively noted for their high catalytic activity towards hydrocarbons for CO 2 RR [26][27][28][29] . However, selective production of CH 4 using Cu catalysts can be difficult owing to the sluggish eight electron transfer steps of CO 2 -to-CH 4 and the inevitable C-C coupling process on bulk Cu 30,31 .
Single-atom catalysts (SACs) with adjustable and isolated active sites have been applied to electrochemical catalysis [32][33][34] , showing potential in expelling C-C coupling in CO 2 RR. Particularly, the unique 3d transition-metal-four nitrogen (M-N 4 ) configuration consisting of an isolated single metal atom coordinated by four N atoms in carbon matrix has been demonstrated to be favorable for CO 2 electroreduction 13,25,35 . However, the performance of Cu-N 4 materials in methane production is still unsatisfactory, specifically, manifesting selectivity to CO at a less cathodic potential (<−1.0 V versus reversible hydrogen electrode, vs. RHE) [36][37][38] and sluggish kinetics to CH 4 at a more cathodic potential 39 . Previous theoretical calculations revealed that relatively weak adsorption of CO* over Cu-N 4 results in feasible CO* desorption instead of further protonation to CH 4 40 . As such, modulation of the electronic structure of Cu-N 4 sites, by substituting relatively weak electronegativity functionalities, e.g., boron (B), is promising to facilitate the adsorption of key intermediates and thus steer the CO 2 RR selectivity towards CH 4 .
Here, we firstly conducted theoretical simulations to predict the adsorption of intermediates and thermodynamic trend over a series of sites with different concentration of boron in Cu-N x B y . It showed that, by introducing boron atoms into the coordination structure of Cu-N 4 , the binding of CO* and CHO* intermediates are promoted significantly, indicating the facile generation of CH 4 . Inspired by the theoretical predictions, we then managed to manipulate the nearest neighbor structure of isolated Cu sites with boron dopant (BNC-Cu). Comprehensive analysis of XAS data revealed that atomically dispersed Cu atoms in BNC-Cu possessed the B-doped Cu-N x structure (Cu-N x B y ) mainly in the form of Cu-N 2 B 2 . Electrochemical measurements displayed a huge boost in CH 4 production, with BNC-Cu showing a high CH 4 Faradaic efficiency (FE) of 73% and partial CH 4 current density (j CH4 ) of −292 mA cm − 2 at −1.46 V vs. RHE. The divergence in CO 2 RR performance between boron-doped and undoped Cu sites validated the theoretical predictions, demonstrating that doping B into Cu-N 4 structure serves as an effective way to enhance deep reduction activity in CO 2 RR. We further studied several Cu-N 2 B 2 -containing sites with the application of twodimensional (quasi) activity and selectivity map along with reaction phase diagram, and the microkinetic simulation of selectivity also matched well with experimental data.

Results and discussion
Thermodynamic trend prediction of CO 2 -to-CH 4 on Cu-N x B y We first built a series of Cu-N x B y structures with different boron concentrations ( Supplementary Fig. 1) and performed density functional theory (DFT) calculations to obtain the adsorption energies of all intermediates in CO 2 -to-CH 4 ( Supplementary Fig. 2). Based on the scaling relations shown in Supplementary Fig. 3, the adsorption free energy of intermediate CHO*, G ad (CHO*), was chosen as the primary descriptor to evaluate the (thermodynamic) activity trend of CO 2 RR to CH 4 . The trend was established considering optimal paths with globally minimal (limiting) energies.
As schematically shown in Fig. 1a, a reaction path consists of many elementary steps. For given paths (black and red), they are not limited by the steps at the beginning of reaction. Instead, r A and r C are the limiting factors. All competing reaction paths were enumerated by the CatRPD code 41 , prior to energetic comparison, as explicitly described in our previous work 41,42 . Briefly, for the internal comparison of a path, the most difficult steps were first determined, such as the r A , r B , and r C from these three different paths in Fig. 1a. Then, the external comparison between all paths can derive the optimal step and energy, defined as ΔG RPD -limiting step and energy. In this case, blue path consisting of r B with the global minimum ΔG is considered as the actual path while r B is ΔG RPD -limiting step. As the ΔG RPD -limiting steps and energies are evolutional with the change of descriptor value, the (quasi) activity and selectivity trend for CO 2 RR can be established, named as reaction phase diagram (RPD) analysis.
Compared to pristine Cu-N 4 , the boron dopants can effectively enhance the reactivity of Cu sites and the adsorption energies. For instance, the more stable COOH* and CHO* can promote the protonation of CO 2 and CO, resulting in promoted production of CH 4 over Cu-N x B y sites at various applied potential ( Fig. 1b and Supplementary  Fig. 4). Projected density of states (PDOS) of CHO* adsorbed on different Cu-N x B y sites explain well more stable adsorption energies ( Fig. 1c-g and Supplementary Fig. 5). Different from the localized 2p-C states of adsorbed CHO* on Cu-N 4 (Fig. 1c), it is more delocalized as the C atom (of CHO*) simultaneously binds with Cu and B at Cu-N x B y sites, which indicates much stronger electronic interaction between CHO* and Cu-N x B y sites. Besides, the hybridized peaks (i-x, as marked in Fig. 1d-g and Supplementary Fig. 5, respectively) below Fermi level show the strong electronic resonance between 2p-C (of CHO*) and 3d-Cu/2p-B of Cu-N 3 B, Cu-N 2 B 2 , Cu-NB 3 , and Cu-N 4 B 4 , suggesting enhanced reactivity on different Cu-N x B y sites. This is also supported by analyzing the crystal orbital Hamilton population (COHP) (Supplementary Fig. 6). The COHP analysis showed that, the intermediate CHO* adsorbs on CuN 4 through C atom bonding with Cu (C-Cu), where the antibonding states of C-Cu bond is partially occupied. With B substituting in Cu-N 4 motif, the C-B bonds are much stronger than C-Cu, indicating the enhancement of CHO* adsorption at Cu-N x B y sites. Besides, the marked peaks below Fermi level (b-e) refer to the strong electronic resonance between 2p-C (of CHO*) and 2p-B of Cu-N x B y sites. Both PDOS and COHP analysis indicated that B substituting in Cu-N 4 motif can promote the adsorption of key intermediates, providing an opportunity to boost the intrinsic activity of CO 2 -to-CH 4 .

Synthesis and characterization of Cu single atom catalysts
Encouraged by these promising theoretical findings, we sought to construct an isolated Cu catalyst in which boron and nitrogen atoms are coordinated with Cu center. The lack of molecules containing Cu-N and Cu-B bond simultaneously limited the ball milling synthesis of a well-defined catalyst with a precise coordination configuration 43,44 . Hence, we managed to introduce B via a one-pot carbonization process (see details in "Methods"). Nevertheless, the high-temperature pyrolysis process may lead to the co-existence of several structures containing Cu-N x B y motifs. We also prepared the Cu-N 4 catalyst (NC-Cu) similarly but without boron precursor (Supplementary Figs. 7 and 8). Transmission electron microscopy (TEM) and scanning electron microscopy (SEM) images indicated that the BNC-Cu catalyst possessed a tubular structure with a mean external diameter of~100 nm (Fig. 2a, b). X-ray diffraction (XRD) pattern and Raman spectrum only exhibited fingerprint peaks corresponding to graphitic carbon without Cu/CuO x /CuB x signals, suggesting the formation of highly dispersed Cu species (Supplementary Figs. 9 and 10). The uniform distribution of Cu, B, C, and N was verified by energy-dispersive X-ray spectroscopy (EDS) elemental mapping ( Supplementary Fig. 11), confirming the formation of B and N co-doped carbon support. The bright spots marked with red circles in high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) image represented isolated Cu atoms dispersed in BNC support (Fig. 2c). Cu particles were not observed in HAADF-STEM images ( Fig. 2c and Supplementary Fig. 11). Inductively coupled plasma atomic emission spectroscopy (ICP-AES) analysis revealed Cu loadings to be 2.5 wt.% for BNC-Cu and 2.7 wt.% for NC-Cu, respectively. X-ray photoelectron spectroscopy (XPS) results indicated the elemental composition of BNC-Cu with N: 17.2%, C: 38.0%, B: 16.9% (Supplementary Fig. 12a), demonstrating heavy doping of N and B into carbon matrix with similar doping-level. The high-resolution XPS spectrum for B 1 s ( Supplementary Fig. 12b) suggested the existence of B-C, B-N, and B-O structure. N 1 s spectrum ( Supplementary Fig. 12c) also displayed obvious component contribution corresponding to N-B [45][46][47] . The detailed information for high-resolution XPS spectra of NC-Cu could be found in Supplementary Fig. 13. Comparison between N 1 s spectra of BNC-Cu and NC-Cu ( Supplementary Fig. 14) also indicated that different N species were generated from general interaction between B and N atoms.
We further employed X-ray absorption near-edge structure (XANES) spectroscopy to accurately identify the nitrogen and boron species in BNC-Cu. The B, N, and C K-edges XANES (Fig. 2d, e and Supplementary Fig. 15) for BNC-Cu exhibited the fingerprint peaks from B-C, B-N, and C-N, demonstrating the strong orbital interaction, via direct bonding, of arbitrary pairs from B, N, and C atoms 45,48,49 . Fourier transform infrared (FTIR) spectrum of BNC-Cu ( Supplementary Fig. 16) also confirmed the existence of similar structures B-N, C-N, and B-C 50 . Meanwhile, the broadened N K-edge spectrum was associated with multiple forms of N species, corresponding well with XPS results. These results strongly supported that boron and nitrogen dopants were integrated into carbon matrix simultaneously, and interaction generally existed between such co-doped boron and nitrogen atoms. We then conducted XAS measurements (Fig. 2f-h

and Supplementary
Figs. [17][18][19][20] of Cu K-edge to reveal the geometric/electronic structure of isolated Cu center in BNC-Cu and NC-Cu. The Fourier transform extended X-ray absorption fine structure (FT-EXAFS) (Fig. 2g) of both samples showed only one main peak with similar radial distance at 1.46 Å which was ascribed to Cu-N/Cu-B scattering path, different from typical Cu-O coordination at 1.57 Å in CuO. No characteristic peak of Cu-Cu contribution at 2.20 Å could be seen, confirming the atomic dispersion of Cu atoms on BNC support. To further investigate the first layer atomic coordination environment of Cu single atoms on BNC and NC support, we performed wavelet transform (WT) of the Cu K-edge EXAFS oscillations ( Fig. 2f and Supplementary Fig. 19). The WT spectra for Cu foil and CuO showed intensity peaks at 7.90 and 6.40 Å -1 , attributed to typical Cu-Cu and Cu-O coordination, respectively. However, for NC-Cu catalyst, the WT spectrum showed one dominated peak, with smaller k value than those of Cu-Cu and Cu-O coordination, indicating the sole formation of Cu-N bonds (Fig. 2f). We further observed the slight WT peak shift of BNC-Cu sample compared to NC-Cu, as shown in Fig. 2f, suggesting the subtle change of atomic coordination environment between BNC-Cu and NC-Cu. A positive shift of absorption edge in BNC-Cu compared with NC-Cu, as shown in Supplementary  Fig. 17, indicated that the electronic structure of Cu atoms changed with the doping of B atoms. EELS point spectrum obtained at the edge of substrate (white box in Supplementary Fig. 21) on a small area, containing three Cu single atoms (bright dots), showed the colocation of B and N atoms around single Cu atoms within the carbon matrix. Considering the angstrom resolution of the electron probe, the signals in EELS point spectrum comes from the closest neighboring atoms of Cu single atoms 51  with NC-Cu ( Supplementary Fig. 13), due to the general interaction between B and N, also indicated the ubiquitous nearby boron atoms around nitrogen. Thus, the more positive Cu valence in BNC-Cu was very likely due to boron atoms that substituted nitrogen atoms and directly coordinated with Cu. Hence, taken above together, we could speculate that boron has been incorporated into the first coordination shell of isolated Cu, forming Cu-N x B y motif. Furthermore, we fitted the EXAFS spectra of BNC-Cu and NC-Cu to illustrate the most possible local configuration of isolated Cu center, as shown in Fig was used as electrolyte (see details in "Methods"). We performed bulk electrolysis over BNC-Cu under different current density with potential range of -0.93 to −2.06 V vs. RHE (Fig. 3a) and over NC-Cu with potential range of −1.18 to −1.96 V vs. RHE (Supplementary Fig. 25).
Gaseous and aqueous products were quantified by gas chromatography (GC) and 1     similar applied potentials. In addition, B-doped Cu-N x active site in BNC-Cu demonstrated significant lower overpotential for CO 2 -to-CH 4 conversion compared to Cu-N 4 by in situ differential electrochemical mass spectrometry (DEMS) analysis ( Fig. 3b and Supplementary  Fig. 28). More specifically, defining the potential, where S/N of m/ z = 15 signal is 5, as the onset potential for CH 4 production, we could extrapolate the CH 4 onset potential of −1.11 V vs. RHE for NC-Cu whereas that of BNC-Cu was −0.82 V vs. RHE (Supplementary Fig. 29). Such difference indicated that introduction of boron ligand into Cu-N 4 motif did further promote its intrinsic activity for CO 2 -to-CH 4 . For BNC-Cu, the hydrogen evolution reaction (HER) firstly dominated the cathodic reaction under low applied potentials, and then retreated dramatically with increased bias while CO 2 RR became the main cathodic reaction (Fig. 3c). The FE CH4 of BNC-Cu kept above 60% at a wide potential range from −1.23 to −1.83 V vs. RHE. The highest FE CH4 of 73% reached at −1.46 V vs. RHE with a j CH4 of −292 mA cm -2 , as shown in Fig. 3d. The j CH4 continuously grew with increasing cathodic potential, reaching to an impressive value of −462 mA cm -2 at −1.94 V vs. RHE. In stark contrast, NC-Cu showed very mediocre performance for CH 4 production. The FE for CH 4 rarely surpassed 30% at tested potential range while the FE for H 2 basically kept above 60% (Supplementary Fig. 25 (Fig. 3e). Post-analysis of the cycled catalyst further demonstrated the structural stability of the formed Cu-N x B y geometry for CO 2 -to-CH 4 (Supplementary Fig. 20 and Supplementary Table 1). For the concern of intensely dynamic aggregation of atomic-dispersed Cu species to clusters during electrolysis, as proposed and demonstrated by other researchers 56,57 , in situ XAS experiments for BNC-Cu and NC-Cu at Cu K-edge under different applied potentials were also carried out. No visible signal of Cu-Cu coordination could be seen in both catalysts when applied with different negative potentials ranging from −0.3 to −1.5 V vs. RHE ( Fig. 3f and Supplementary Fig. 35), indicating that the atomically dispersed Cu atoms remained as dominant Cu species during electrolysis. The above results revealed that BNC-Cu catalyst performed much better than the boron-free counterpart NC-Cu for converting CO 2 to CH 4, validating the theoretical prediction results.

Insight into CO 2 -to-CH 4 conversion mechanism
To further elucidate the CO 2 RR to CH 4 mechanism over BNC-Cu, in situ attenuated total reflectance surface-enhanced infrared absorption spectroscopy (ATR-SEIRAS) measurements were conducted for BNC-Cu, as well as the NC-Cu control ( Fig. 4a and Supplementary Fig. 36). As depicted in Fig. 4a, on applying potential from −0.8 to −1.5 V vs. RHE over BNC-Cu, two peaks could be well noted. Both peaks showed a slight red-shift with increasing applied potential, indicating that such observed signals originated from the in situ formed CO 2 RR intermediates 58 . We ascribed the detected two peaks at around 2100 and 1720 cm -1 to CO* 59,60 and the key intermediate to produce CH 4 [21][22][23]61 , CHO* 62 , respectively. On scanning to a negative potential, with the gradual disappearance of CO* signal, we observed the accumulation of CHO* intermediate, revealing the rapid transformation of CO* into CHO* at high applied overpotential. On the other hand, no visible CHO* peak was observed for NC-Cu ( Supplementary Fig. 36), while only the CO* fingerprint peak appeared and the intensity remained almost constant over the tested potential range, revealing higher barrier for CO* hydrogenation to CHO* on NC-Cu. Such a difference on CO*-to-CHO* transition may explain the CO 2 RR performance gap between BNC-Cu and NC-Cu. Besides, the ATR-SEIRAS displayed another two strong broad peaks at around 1650 and 3400 cm -1 , which could be ascribed to the interfacial H 2 O 63,64 . We noted the intensity of these peaks increased with increasing applied potential, suggesting that water molecules tended to adsorb on the surface of BNC-Cu during CO 2 RR (Fig. 3a) 65 . Such readily adsorbed water molecules could provide adequate hydrogen source for CH 4 production, but this trend was not obvious for NC-Cu (Supplementary Fig. 36).
The observed signals of CHO* further corroborated the use of CHO* as the target in the aforementioned computational screening over a series of Cu-N x B y sites ( Supplementary Figs. 2 and 3). Based on the resolving results of Cu center in BNC-Cu, it is reasonable to speculate the boost of CH 4 production was very likely due to the dominant Cu-N 2 B 2 structure, and we thus performed further theoretical investigations to understand such high selectivity towards CH 4 with barrier calculations in consideration of several structures containing Cu-N 2 B 2 motif (Cu-N 2 B 2 −1, Cu-N 2 B 2 −2, Cu-N 2 B 2 −3, Cu-N 4 B 4 −1, and Cu-N 4 B 4 −3, illustrated in Supplementary Figs. 2 and 3). On the basis of one dimensional (quasi) activity map using G ad (CHO*) as the primary descriptor (Fig. 1b), we chose G ad (CO*) as the second descriptor as the adsorption strength of CO* is tightly correlated with CO production, thus forming a two-dimensional (quasi) activity and selectivity map that clearly shows the selectivity trend of Cu-N 4 , Cu-N 2 B 2 , and Cu-N 4 B 4 sites (Fig. 4b). At −1.46 V vs. RHE, the maximum FE CH4 of BNC-Cu observed in experiments, HER is thermodynamically favorable at sites with weak surface reactivity, for instance, Cu-N 4 site. It again demonstrates how corporation of B atoms into Cu-N 4 structure augments CO 2 RR to CH 4 . However, CO 2 RR is the main reaction on sites with strong adsorption of intermediates, where relatively weak CO* adsorption is prone to desorb and produce CO, otherwise to CH 4 production. Accordingly, the site Cu-N 2 B 2 −2, as well as Cu-N 4 B 4 −3, are more likely to contribute in CO 2 to CH 4 process among other structures.
To further make comparison between Cu-N 2 B 2 −2 and Cu-N 4 B 4 −3, electrochemical barriers of CO 2 RR and HER over the two structures were calculated and shown in Fig. 4c-f. The electrochemical barriers of proton-coupled electron transfer reactions were calculated via a capacitor model 66,67 . According to "charge-extrapolation" method 67,68 , the amount of electron transfer (Δq) from water to the electrode surface is linearly correlated with the relative work function (Φ) at the initial states (IS), transition states (TS), and final states (FS). Figure 4c-e show such linear correlations of the critical steps for products selectivity. Illustrated by Fig. 4c, the TS of CO* protonation step to CHO* is IS-like, giving rise to a small charge transfer coefficient (β). A median TS of the ΔG RPD -limiting step for CO 2 RR, CH 3 * protonation (Fig. 4d), leads to a moderate β. However, the β of Heyrovsky step, the ΔG RPDlimiting step for HER, is larger due to its FS-like TS (Fig. 4e). It accurately predicts the high FE of H 2 at very negative potentials. Based on calculated kinetic barriers, a detailed free energy landscape of CO 2 RR over Cu-N 2 B 2 −2 and Cu-N 4 B 4 −3 sites is shown in Fig. 4f (HER in subfigure). The CO* protonation has a barrier of 0.43 eV on Cu-N 2 B 2 −2, while such process is barrier-less on Cu-N 4 B 4 −3. However, for CH 3 * protonation, the larger barrier on Cu-N 4 B 4 −3 of 1.40 eV than that on Cu-N 2 B 2 −2 of 0.48 eV supports well the significance of global energy optimization. Specifically, Cu-N 2 B 2 −2 shows lower barriers of the most energetically difficult steps than Cu-N 4 B 4 −3 for CO 2 RR, which is consistent with the (quasi) activity map in Fig. 4b. Thus, Cu-N 2 B 2 −2 site is expected to show prominent CH 4 production during CO 2 RR. The following microkinetic modeling was conducted for Cu-N 2 B 2 −2 at −1.46 V vs. RHE (Fig. 4g). The similarity of selectivity between computational and experimental data inversely suggested that Cu-N 2 B 2 −2 sites are probably the dominant sites among various possible Cu-containing sites. The experimental FE CH4 value was 75%, which was lower than the simulated result of 80%. This discrepancy was due to the existence of other sites that were present in small proportions and did not have as high a methane productivity. The comparable product distribution also inversely validated our speculation that the majority of Cucontaining sites was Cu-N 2 B 2 , resolved by multiple XAS-relative studies. Compared with the pristine Cu-N 4 , the enhanced surface reactivity and appropriate adsorption energies over Cu-N 2 B 2 −2 is thus proposed as the chemical origin of the high selectivity towards CH 4 for as-synthesized BNC-Cu catalyst.
In summary, by virtue of thermodynamic trend and global energy optimization analysis, refined manipulation of experimental synthesis and characterizations, theoretical microkinetic modeling, this work showcased a general way to design Cu single atom catalysts towards methane production via rationally regulating the nearest coordination environment. By modifying the Cu-N 4 sites with partial B substitution, the enhanced adsorption of CO* and CHO* intermediates are theoretically proven to be beneficial for CH 4 generation. Experimentally, asobtained B-doped Cu-N x active sites exhibited a high CH 4 FE of 73% at −1.46 V vs. RHE and a maximal j CH4 of −462 mA cm -2 at −1.94 V vs. RHE, respectively. This work implicates a valuable avenue for bolstering selectivity of Cu towards a specific product.

Computational details
All density functional theory (DFT) calculations in this work were performed via the Vienna Ab-initio Simulation Package (VASP) 69,70 . The generalized gradient approximation (GGA) 71 with the revised Perdew-Burke-Ernzerhof (rPBE) functional 72 (GGA-rPBE) was employed to describe the electron interactions. The projected augmented wave (PAW) 73,74 was employed to describe the valence electrons with a plane wave basis sets and the kinetic energy cutoff was set to 400 eV. Structural optimizations were performed with the residual force and electronic energy differences smaller than −0.05 eV Å -1 and 10 -5 eV, respectively. In addition, dispersion effects have been taken into account by DFT-D3 method. Locating transition states were conducted by the climbing image nudged elastic band (CI-NEB) 75 and dimer 76,77 methods, where the convergence force was set as smaller than 0.1 eV Å −1 . Moreover, a Gaussian smearing with a width of 0.2 eV and a Monkhorst-Pack k-point mesh grid of 2 × 2 × 1 were used. For atomic model construction, Cu-N 4 was embedded in a graphene structure of 17.16 × 14.87 Å 2 as the active center of NC-Cu (Supplementary Fig. 1), where a vacuum of 15 Å was introduced to avoid the interaction of adjacent layers. Several typical structures with different amounts of boron incorporated were also optimized (Supplementary Fig. 1). The adsorption energy (E ad ) of intermediates refers to the energies of gas-phase CO(g), H 2 (g), and H 2 O(g), the chemical potential of single H, O, C atom is The intermediates adsorption energy (E ad ) can be computed by the following equation: where E bare and E tot are the energies of bare catalysts and the catalysts with adsorbates, respectively. The numbers of carbon, hydrogen and oxygen atoms in intermediates are described by coefficients α, β, and γ. Furthermore, free energy corrections were conducted in this work at the temperature of 298 K, with the scheme described in the previous work 78 .
The formation energy of CuN x B 4-x sites (x = 0-3) are computed via the following equation: where n i and μ i refer to the number and chemical potential of element i, respectively. The energies of Cu and C are referenced to the Cu bulk and graphene, respectively. Besides, the chemical potentials of N and B are calculated via the energies of HNO 3 , H 3 BO 3 , NH 2 CN, and H 2 O.

The elementary steps of electrochemical CO 2 RR and HER
We have considered three whole reactions: In which several elementary reactions may be involved, as listed in the following: COOH * + H 2 O ! HOCOH * + ðOH À À e À Þ ðR2Þ COOH * ! CO * + ðOH À À e À Þ ðR4Þ Theoretical XANES spectrum calculations The Cu K-edge XANES simulation was conducted with the FDMNES code in multiple scattering mode (Green) using the muffin-tin potential. The energy dependent exchange-correlation potential was calculated in the real Hedin-Lundqvist scheme, and then the spectra convoluted using a Lorentzian function with an energy-dependent width to account for the broadening due both to the core-hole width and to the final state width 79,80 . Parameter optimization was performed by comparing the theoretical and experimental spectra to acquire the most appropriate convolution parameters. The calculated models were built based on DFT calculations to avoid manual bias.

Chemicals
All chemicals were used as received without further purification. The NC-Cu catalyst was synthesized similarly with 0.1 g copper (II) nitrite trihydrate input and without boric acid. The pure-BNC substrate was also synthesized without copper precursor, while pure NC substrate was synthesized without both coppern precursor and boric acid.

Characterization techniques
The as-synthesized BNC-Cu and NC-Cu were characterized by various analytical techniques. X-ray diffraction (XRD) was performed on a Philips X'Pert Pro Super diffractometer with Cu-Kα radiation (λ = 1.54178 Å). The morphology of the samples was observed by scanning electron microscopy (SEM, Zersss Supra 40) and transmission electron microscopy (TEM, Hitachi H-7650). HAADF-STEM images, energy-dispersive X-ray Spectroscopy (EDS) elemental mapping, and electron energy loss spectroscopy (EELS) were carried out on JEOL ARM-200F field-emission transmission electron microscope operating at an accelerating voltage of 200 kV using Mo-based TEM grids. Raman spectra were taken on a Raman microscope (Renishaw®) excited with a 785 nm excitation laser. X-ray photoelectron spectroscopy (XPS) measurements were performed on a VG ESCALAB MK II X-ray photoelectron spectrometer with Mg Kα = 1253.6 eV as the exciting source. Soft X-ray absorption spectra (B K-edge, N K-edge, and C K-edge) were carried out at the Catalysis and Surface Science Endstation at the BL11U beamline in the National Synchrotron Radiation Laboratory (NSRL) in Hefei, China.

Electrochemical measurements
The electrochemical measurements were conducted with an electrochemical workstation (CHI 660E, Shanghai CH Instruments). The Ag/ AgCl wire in saturated KCl solution was adopted as the reference electrode, and the counter anodic reaction was oxygen evolution reaction over a Ni foam. All potentials were converted to the RHE reference scale using the relation E RHE = E Ag/AgCl + 0.197 + pH × 0.059. Solution resistance was determined by potentiostatic electrochemical impedance spectroscopy at frequencies ranging from 0.1 Hz to 200 kHz and compensated by 85%. To prepare a catalyst cathode in a flow reactor, the catalyst ink was prepared at first with a constant composition ratio for each sample, 10 mg of catalyst, mixed with 25 μL of ionomer in 2 mL of isopropanol. Such ink was then airbrushed onto a 2 cm × 1.5 cm carbon gas diffusion layer (39BC) under the heating of a heat stage set in 80°C. The mass loading for each catalyst was determined by weighing the mass of carbon paper before and after sprayed with catalyst ink and controlled to be~1.0 mg cm −2 . Such two electrodes were then placed on opposite sides of two polytetrafluoroethylene (PTFE) sheets with 0.4 cm × 1.5 cm channels. The geometric surface area of catalysts was controlled as 0.6 cm 2 . A Nafion 115 membrane (Fuel Cell Store) was sandwiched between the two PTFE sheets to separate the chambers. On the cathode side, a titanium gas flow chamber supplied 30 s.c.c.m. CO 2 (monitored by an Alicat Scientific mass flow controller). 0.5 M KHCO 3 electrolyte was pumped through cathode chambers with a constant rate of 0.75 mL min −1 , while 1 M KOH served as anodic electrolyte was circulated around the anode with a rate of 25 mL min −1 .

CO 2 reduction product analysis
To quantify the gas products obtained during CO 2 electrolysis, pure CO 2 gas was delivered into the cathodic compartment at a constant rate and vented into a gas chromatograph (PerkinElmer Clarus® 690) equipped with a thermal conductivity detector and a flame ionization detector. The liquid products were quantified using a 400 MHz NMR spectrometer. Typically, after electrolysis, 600 μL electrolyte was mixed with 100 μL D 2 O (Sigma Aldrich, 99.9 at.% D) and 0.05 μL dimethylsulfoxide (Sigma Aldrich, 99.9%) as internal standard.

Ex situ and in situ XAS experiments
The ex situ X-ray absorption spectroscopy (XAS) spectra of Cu K-edge were obtained using beamline 44 A of Taiwan Photon Source (TPS) at National Synchrotron Radiation Research Center, Taiwan. All XAS data of Cu K-edge were collected in fluorescence mode using 7-element SDD detector and the incident photon energy were calibrated using standard Cu foil. In situ XAS spectra of Cu K-edge were obtained using beamline BL11B and BL14W at the Shanghai Synchrotron Radiation Facility (SSRF), Shanghai advanced Research Institute, Chinese Academy of Sciences. All XAS data of Cu K-edge were collected in fluorescence mode and the incident photon energy were calibrated using standard Cu foil. A self-design organic glass electrochemical cell was set in a three-electrode configuration and employed for our in situ XAS experiments. A rectangular orgonic glass cap was used to cover the cell and to keep the cell at fixed position of optical path. Several holes on the cap are used for CO 2 bubbling and ensuring a fixed distance between working and reference electrodes for all experiment. A graphite rod and an Ag/AgCl electrode were used as the counter electrode and reference electrode, respectively. The working cell has flat walls with a single circular hole of 2 cm in diameter as a window of contact between electrolyte and catalysts, and a beam of synchrotron radiation X-ray light irradiated within the circular area during the in situ XAS experiments. To prepare a catalyst working electrode, 18 mg of catalyst, mixed with 45 μL of ionomer in 4 mL of isopropanol, such ink was then airbrushed onto a 2.5 cm × 2.5 cm carbon gas diffusion layer (39BC) with a mass loading of~1.0 mg cm −2 . Catalyst coated carbon paper was in contact with a slip of copper with the catalyst layer facing inward. Then 20 mL 0.5 M KHCO 3 solution pre-saturated with CO 2 was poured into the cell. The solution was not stirred and CO 2 was bubbled into the solution bottom through the hole on that orgonic glass cap during the experiment. The flow rate of CO 2 was 10 s.c.c.m. monitored by an Alicat Scientific mass flow controller. The cell was connected to an electrochemical station by making electrical contact to the copper tape slip that protruded from the side of the working cell. Before the in situ XAS experiments, XAFS spectra were recorded at different positions on the electrode to check the homogeneity of the catalyst. During the in situ XAS experiments, the potential on working electrode started from ocp to a series of cathodic potentials, and back to ocp. At each potential, the system was allowed to equilibrate for 10 min before recording a spectrum, then scans at the Cu K-edge were recorded. Data reduction, data analysis, and EXAFS fitting for XAS analysis in this work were performed with the Athena, Artemis, and IFEFFIT software packages. For quantitative analysis, phase shifts and backscattering amplitudes were generated by the FEFF calculations based on crystal structures of Cu, and were then calibrated through performing the FEFFIT of the EXAFS data of the reference samples, mainly to obtain the amplitude reduction factor (S 0 2 ) values. With S 0 2 known, the EXAFS data of the catalyst materials were fitted with such generated phase shifts and amplitudes.

In situ DEMS experiments
The "probe-type" DEMS was applied for the detection of volatile CH 4 and H 2 produced during the CO 2 RR, as well as the reactant CO 2 consumed. The whole tests were conducted in a flow-cell system with a gas diffusion electrode as working electrode. To prepare a catalyst working electrode, 10 mg of catalyst, mixed with 25 μL of ionomer in 2 mL of isopropanol, such ink was then airbrushed onto a 2 cm × 1.5 cm carbon gas diffusion layer (39BC) with a mass loading of 1.0 mg cm −2 . Catalyst coated carbon paper was in contact with a slip of copper with the catalyst layer facing inward to the electrolyte. The effective area of the working electrode is a circle of 1 cm in diameter. A sampling probe approached the working electrode at a distance of ca. 20 μm, and a peristaltic pump replaced the solution near the working electrode at a flow rate of 1.25 mL min -1 . The onset potential for each product is defined as the potential where the S/N of corresponding m/z signal is 5.
In situ attenuated total reflection surface-enhanced IR absorption spectroscopy (ATR-SEIRAS) measurements In situ ATR-SEIRAS spectrum was gathered by a FT-IR spectrometer (Thermo Scientific Nicolet iS50) equipped with MCT-A detector. The catalyst inks were prepared by mixing 10 mg electrocatalysts, 5 mL ethanol, and 25 μL of ionomer. 10 μL of ink solution was dropped onto the central area (confined by an O-ring with Φ = 8 mm) of an Au film deposited on the basal plane of a hemicylindrical Si prism by evaporation. The Si prism was assembled in a spectro-electrochemical cell with Pt wire as a counter electrode, Ag/AgCl wire in saturated KCl solution as reference electrode, and 0.5 M KHCO 3 solution presaturated and continuously bubbled with CO 2 as electrolyte. All spectra are collected at a resolution of 4 cm −1 and each single-beam spectrum is an average of 200 scans. A CHI 660e electrochemistry workstation (Shanghai CH Instruments, Inc.) was used for potential control.